# Dissolve the district polygons to form new polygon of Bihar and EUPlibrary(sf)India_aoi_sf_dis=st_union(India_aoi_sf)
3 Causal Random Forest Model to Get Individual Treatment Effects
4 Point-geocoded data
4.1 Gridded data input variables
The first strategy is to translate all variables to the grid. This involves interpolation across space and using new variable names. In this case, instead of gender being a dummy, you use a proportion of female or male after interpolation.
4.3.2 Spatial Bayesian Geostatistical Gaussian Process Model
If one is interested in calculating other measures other than the predicted value (for example, the probability of exceeding some amount), then a Bayesian gaussian process model is the best alternative in that using Markov Chain Monte Carlo simulations one can use a probabilistic assessment.
### Bayesian kriging LDS_small <- LDS[sample(1:nrow(LDS),800),] coords=dplyr::select(LDS_small,O.largestPlotGPS.Longitude,O.largestPlotGPS.Latitude)coords=as.matrix(coords)# The public version of the data has duplicated coordinates# We need to jitter these because spatial Bayesian kriging requires unique coordinates.library(geoR)coords=jitterDupCoords(coords,min=2,max=10)coords=as.matrix(coords)# Bayesian models take much time to render. We sample 1000 observations to showcase the approachlibrary(spBayes)n.samples=1000t1 <-Sys.time()r <-1n.ltr <- r*(r+1)/2priors <-list("phi.Unif"=list(rep(1,r), rep(10,r)), "K.IW"=list(r, diag(rep(1,r))), "tau.sq.IG"=c(2, 1))starting <-list("phi"=rep(3/0.5,r), "A"=rep(1,n.ltr), "tau.sq"=1)tuning <-list("phi"=rep(0.1,r), "A"=rep(0.01, n.ltr), "tau.sq"=0.01)cf.sowing.sp <- spBayes::spSVC(yield_kgperha~1, data=LDS_small,coords=coords,starting= starting,tuning=tuning,priors=priors,cov.model="exponential",n.samples=n.samples,n.omp.threads=15,svc.cols=c("(Intercept)"))
----------------------------------------
General model description
----------------------------------------
Model fit with 800 observations.
Number of covariates 1.
Number of space varying covariates 1.
Using the exponential spatial correlation model.
Number of MCMC samples 1000.
Priors and hyperpriors:
beta flat.
K IW hyperpriors:
df: 1.00000
S:
1.000
phi Unif lower bound hyperpriors: 1.000
phi Unif upper bound hyperpriors: 10.000
tau.sq IG hyperpriors shape=2.00000 and scale=1.00000
Source compiled with OpenMP, posterior sampling is using 15 thread(s).
-------------------------------------------------
Sampling
-------------------------------------------------
Sampled: 100 of 1000, 10.00%
Report interval Metrop. Acceptance rate: 64.00%
Overall Metrop. Acceptance rate: 64.00%
-------------------------------------------------
Sampled: 200 of 1000, 20.00%
Report interval Metrop. Acceptance rate: 42.00%
Overall Metrop. Acceptance rate: 53.00%
-------------------------------------------------
Sampled: 300 of 1000, 30.00%
Report interval Metrop. Acceptance rate: 22.00%
Overall Metrop. Acceptance rate: 42.67%
-------------------------------------------------
Sampled: 400 of 1000, 40.00%
Report interval Metrop. Acceptance rate: 24.00%
Overall Metrop. Acceptance rate: 38.00%
-------------------------------------------------
Sampled: 500 of 1000, 50.00%
Report interval Metrop. Acceptance rate: 30.00%
Overall Metrop. Acceptance rate: 36.40%
-------------------------------------------------
Sampled: 600 of 1000, 60.00%
Report interval Metrop. Acceptance rate: 30.00%
Overall Metrop. Acceptance rate: 35.33%
-------------------------------------------------
Sampled: 700 of 1000, 70.00%
Report interval Metrop. Acceptance rate: 32.00%
Overall Metrop. Acceptance rate: 34.86%
-------------------------------------------------
Sampled: 800 of 1000, 80.00%
Report interval Metrop. Acceptance rate: 24.00%
Overall Metrop. Acceptance rate: 33.50%
-------------------------------------------------
Sampled: 900 of 1000, 90.00%
Report interval Metrop. Acceptance rate: 34.00%
Overall Metrop. Acceptance rate: 33.56%
-------------------------------------------------
Sampled: 1000 of 1000, 100.00%
Report interval Metrop. Acceptance rate: 23.00%
Overall Metrop. Acceptance rate: 32.50%
-------------------------------------------------
Source compiled with OpenMP, posterior sampling is using 1 thread(s).
-------------------------------------------------
Recovering beta and w
-------------------------------------------------
Sampled: 99 of 251, 39.44%
Sampled: 199 of 251, 79.28%
# Kriginglibrary(terra)library(stars)library(raster)library(gstat) # Use gstat's idw routinelibrary(sp) # Used for the spsample functionlibrary(tmap)library(geodata)# India=gadm(country="IND", level=1, path=tempdir())# plot(India)# India_State_Boundary=subset(India,India$NAME_1=="Bihar")# plot(India_State_Boundary)# library(sf)#India_State_Boundary=st_as_sf(India_State_Boundary)#India_State_Boundary_Bihar_sp=as_Spatial(India_State_Boundary_Bihar)library(spBayes)India_aoi_sf_dis_sp=as_Spatial(India_aoi_sf_dis)LDS_small_sp=SpatialPointsDataFrame(cbind(LDS_small$O.largestPlotGPS.Longitude,LDS_small$O.largestPlotGPS.Latitude),data=LDS_small,proj4string=CRS("+proj=longlat +datum=WGS84"))LDS_small_sp@bbox <- India_aoi_sf_dis_sp@bboxgrd <-as.data.frame(spsample(LDS_small_sp, "regular", n=10000))names(grd) <-c("X", "Y")coordinates(grd) <-c("X", "Y")gridded(grd) <-TRUE# Create SpatialPixel objectfullgrid(grd) <-TRUE# Create SpatialGrid objectplot(grd)
Source compiled with OpenMP, posterior sampling is using 15 thread(s).
-------------------------------------------------
Point-wise sampling of predicted w
-------------------------------------------------
Sampled: 99 of 251, 39.44%
Sampled: 199 of 251, 79.28%
# Standard deviationlibrary(rasterVis)pred.sd=pred.grid["pred.sd"]pred.sd=raster(pred.sd)pred.sd=mask(pred.sd,India_aoi_sf_dis_sp)pred.sd_plot=levelplot(pred.sd,par.settings=RdBuTheme(),contour=TRUE)pred.sd_plot
# Probability of 3 tonspred.prob_3tons=pred.grid["pred.prob_3tons"]pred.prob_3tons=raster(pred.prob_3tons)pred.prob_3tons=mask(pred.prob_3tons,India_aoi_sf_dis_sp)pred.prob_3tons_plot=levelplot(pred.prob_3tons,par.settings=RdBuTheme(),contour=TRUE)pred.prob_3tons_plot
# Probability of 5 tonspred.prob_5tons=pred.grid["pred.prob_5tons"]pred.prob_5tons=raster(pred.prob_5tons)pred.prob_5tons=mask(pred.prob_5tons,India_aoi_sf_dis_sp)pred.prob_5tons_plot=levelplot(pred.prob_5tons,par.settings=RdBuTheme(),contour=TRUE)pred.prob_5tons_plot
4.3.3 Spatial Bayesian Geoadditive Model
5 Areal data
5.1 Areal centroid random field gaussian process model
One can get the centroid of each polygon, e.g., then fit a geostatistical model. This can then be gridded across the area of interest.
5.2 Areal-point downscaling
5.2.1 Markov Random Field (MRF) Geoaddive Structured and Unstructured Spatial Model